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I. 



INTRODUCTION 



Understanding of the characteristics of the airflow over 
an airfoil is of paramount importance to the airfoil 
designer. Two methods are currently available which give 
accurate results. The first is the use of wind tunnel tests. 
The drawbacks to this method are cost and time consumption. 
The second is the processing of the Navier-Stokes equations. 
This method's drawbacks are the requirements and expense of 
using supercomputers due to the extensive calculations and 
storage requirements. There is still a need to come up with 
an inexpensive, fast and accurate engineering tool to compute 
airfoil flows. 

Several methods have been derived to accomplish this end. 
But the most promising is the Viscous-Inviscid Interaction 
method. The outer flow is computed using inviscid flow 
equations, and the inner flow is computed using Prandtl ' s 
boundary layer equations. The key to this method is the 
extent of interaction between the inner and outer flows. 

The purpose of this thesis is to evaluate the capability 
of the viscous-inviscid interactive aircode developed by 
Tuncer Cebeci and associates at the Douglas Aircraft Company 
[Ref. 1]. This computer program was applied to four airfoils 
with various angles of attack and Reynolds numbers. The 
computer results were then compared to previously reported 
experimental results. 
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The conservation of mass and momentum are summarized in 



Chapter 2, inviscid flow calculations are discussed in 
Chapter 3, and viscous flow equations are described in 
Chapter 4. Viscous calculations are presented in Chapter 5, 
and the specific interaction methods are shown in Chapter 6. 
Finally, in Chapter 7 computer and experimental results are 
compared for the NACA 66 3 -018, 0010 (Modified) and 4412 
airfoils as well as the Wortmann FX 63-137 airfoil. 
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II. FUNDAMENTAL EQUATIONS 



The conservation of mass and conservation of momentum 
provide the foundation for incompressible flow analysis. 
With these fundamental concepts along with appropriate 
assumptions and approximations, working relations for two- 
dimensional, incompressible flow are obtained. 

A. CONSERVATION OF MASS (CONTINUITY) 

The conservation of mass principle states that mass 
cannot be created nor destroyed. Equating this statement to 
a fixed control volume the net mass flow rate into and out of 
the control volume equals the time rate of change of mass 
within the control volume [Ref. 2:p. A-l]. 

Given a control volume, the mass flow rate through one of 
its surfaces is equal to the product of the fluid density, 
the fluid velocity normal to the surface and the area of the 
surface [Ref. 3:p. 29]. 



d mass 

= V • n s (2.1) 

dt 

In 2-D flow the x-component of the mass flow rate at the 
center of the positive x-face, position dx/2 and side length 
dy, is represented by Taylor series expansion [Ref. 2:p. A-2] 



3 



+ . . 



dy. 



( 2 . 2 ) 



d dx £) 2 dx 1 

pu + — (pu) — + (pu)( — } 2 — 

c)x 2 dx 2 2 2! 



As dx approaches zero all higher order terms disappear 
leaving 



pu 



d dx 

+ — (pu) — 
dx 2 



dy. 



(2.3) 



Similarly, the x-component of the mass flow rate at the 
center of the negative x-face, position -dx/2 and side length 
dy, is represented by 



pu 



d dx 

— ( pu) — 

dx 2 



dy. 



(2.4) 



As illustrated in Figure 2.1 the net mass flow through the 
four sides of the 2-D control surface is 



d , dx 


dy - 


3 , dx 


pu - — ( pu) — 


pu + (pu) 


dx 2 




dx 2 


d dy 




d dy 


pv - — (pv) 


dx - 


pv + — ( pv) — 


dy 2 




dy 2 
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0 , dx 

pu - — ( pu) — 

c)x 2 



a d y 

pv + — (pv) 

dy 2 

f 




pv - — ( pv) 

0y 2 



y 



> X 



a dx 

> pu + — (pu) 

ax 2 



Figure 2.1 Mass Flow Through 2-D Control Surface 
[ Ref . 4 : p . 12] 
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which is equal to the rate of change of the mass within the 
control volume 



c)P , , 

— dxdy 

c)t 

Combining (2.5) with (2.6) and simplifying yields 


(2.6) 


- — (pu)dxdy - — (pu)dxdy = — dxdy 
c)x c)y c)x 

Dividing by dxdy and rearranging yields 


(2.7) 


~ + ~(pu) + — (p U ) = 0 
at ax c)y 


(2.8) 


For steady, incompressible flows the continuity 
becomes 


equation 


c)u c)v 

— + — =0 
c)x c)y 


(2.9) 


or in vector form the continuity equation [Ref. 3:p. 


30 ] is 


V . v = 0 

B. CONSERVATION OF MOMENTUM ( NAVIER-STOKES ) 


(2.10) 



The conservation of momentum, Newton's second law of 
motion, states that the rate of change of the linear momentum 
is equal to the sum of the forces applied [Ref. 2:p. B-1J. 
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d 

SF = — (mV) (2.11) 

dx 

As illustrated in Figure 2.2 the two significant forces which 
act on an element of fluid are surface forces which act on 
the surface only, pressure and shear, and body forces which 
affect the mass of the element, such as gravity. Assuming 
moment equilibrium in an element, x^y = x y>t , the 2-D first 
order Taylor series expansion for normal and shear surface 
forces in the x-direction is 



a dx ^ 

x X + (Xxx) Xxx + " ' ( X X X ) — 

c)x 2 c)x 2 



dx 



dy 



a , dy _ a , dy 

X >cy + ( X X y ) “ Xxy + ( X X V ) 

Oy 2 ay 2 



dx 



= ^— ( x xlt ) dxdy + — (x^v)dxdy. (2.12) 

3x c)y 

The body forces per unit mass are represented by 

Fbody = Xi + Yj + zk (2.13) 

such that the x-component of the body force on an element is 

f x ( bod v > = pdxdyl»X. (2.14) 

Combining equations (2.12) and (2.14) the sum of the forces 
in the x-direction is 
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d dy 

Xyy + — " ( X yy ) — ■ 

dy 2 



d dx 

- — (x^) — 
dx 2 



a, 

3x 



<- 



0 dy 

Xyx “ ( X y x ) 

dy 2 



A 



-> 




V 



0 dy 

Xyx + ( X y x ) 

dy 2 



3 dx 

X x y + ( X x y ) 

dx 2 
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X X X ( X XX ) 
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Xyy "" |Xyy) 
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Figure 2.2 Stresses on a 2-D Control Surface 
[ Ref . 4 : p . 15 ] 
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dxdy . 



(2.15) 



SF X 



pX + 



a 

- — ( t XX ) 

Ox 



a 

+ ( X x y ) 

ay 



The rate of change of the linear momentum in the x-direction 
assuming constant mass is mdu/dt. 



du 


au dx 


du 


dy 


au 


a x 


= 





+ 





+ — via the chain rule. 


and — 


dt 


ax dt 


ay 


dt 


at 


at 


ay 













and — = v the x-direction change in linear momentum for 

at 

particle is 



du 3u Ou 3u 

m — = pdxdy(v — + v — + — ). (2.16) 

dt dx Oy Dt 



Substituting equations (2.15) and (2.16) into the x-component 
of equation (2.11) yields 



au 

pdxdy(u — + 

ax 



au au 

V + ) 

ay at 



a a 

P X 1 ( "C xx ) F ( ^ xy 

ax ay 



dxdy. (2.17) 



Now, in order to have the entire equation as a function 
of velocity the normal and shear stresses must be found in 
terms of velocity. 

By assuming a Newtonian fluid [Ref. 2:p. B-5] the shear 
stress is linearly related to the rate of angular deformation 
with fluid viscosity being the proportionality constant. 
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(2.18) 



Ou c)v 

X>cy — U( + ) 

c)y Ox 

The normal stresses are equal, but opposite in direction 
to the pressure when no shear stresses are involved. With 
shear stress from viscosity it is assumed that the normal 
stresses deviate from -P and that the deviation is 
proportional to both a) the rate of linear strain in the 
direction concerned, and b) the rate of volume deformation. 
Therefore, the normal stress in the x-direction [Ref. l:p. B- 
10] is 



Ou 2 Ou Ov 

T*x = - P + 2u( — ) - — U ( — + — )• (2.19) 

Ox 3 Ox Oy 



Applying the conservation of mass, equation (2.9), equation 
(2.19) simplifies to 



Ou 

Txx = - P + 2p( — ) • (2.20) 

Ox 



Substituting equations (2.18) and (2.20) into (2.17) and 
dividing by dxdy yields 



Ou 


Ou c)u 


a 


3u 


Q Ou 0v 


p(u 


+ V + ) 


= px + — (-P + 


2U( — )) + 


U— (— + ) (2.21) 


Ox 


c)y at 


c)x 


Ox 


0y 0y Ox 


After 


multiplying 


and rearranging 


the right 


hand side becomes 
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dp d 2 U d 2 U <) 2 V 

px - — + 2y + y + y 

dx dx 2 dy 2 dydx 

which is also equal to 



dp 


d 2 u 


d 2 u 


d du 


dv 


px - — + 


n- — 


+ y 


+ u— ( — 
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dx 


dx 2 


dy 2 


dx dx 


dy 



Again applying the conservation of mass, equation (2.9), 
equation (2.21) becomes 



du 


du 


du 


dp d 2 u d 2 u 


p(u — + 


V 


+ r-) 


= px - — + y( + 


dx 


dy 


dt 


dx dx 2 dy 2 



With v = y/p = kinematic viscosity and neglecting the body 
force, X, the two dimensional Navier-Stokes, conservation of 
momentum equation for incompressible and constant viscosity 
flow in the x-direction [Ref. 2:p. 436] is 



du 


du 


du 


l dp 


d 2 u 


d 2 u 


+ 


u — + 


v — 


— + 


v ( 


+ 


dt 


dx 


dy 


p dx 


dx 2 


dy 2 



Similarly, in the y-direction the Navier-Stokes equation is 



dv 


dv 


dv 


l dp d 2 v d 2 v 


+ 


u + 


V— - — 


+ V ( + 


dt 


dx 


dy 


p dy dx 2 dy 2 



(2.24) 
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III. STEADY INVISCID FLOW 



Although real fluids are viscous the major effects of 
viscosity are concentrated in a region, or layer adjacent to 
a body. Therefore, analyses of inviscid flow are useful and 
serve as a good approximation to flow outside the boundary 
layer and wake behind the body. 

The justification for applying the results of perfect 
fluid analyses to viscous flows was postulated by Ludwig 
Prandtl in 1904 [Ref. 3:p. 299]. He stated that the effects 
of viscosity on the flow around streamlined bodies at high 
Reynolds numbers are effectively limited to a "thin" boundary 
layer. The characteristic length to judge thinness is the 
distance from the forward stagnation point to the point of 
consideration . 

A. VELOCITY POTENTIAL 

For flow outside the boundary layer it is a great 
advantage to simplify equations and develop a single 
governing equation. With the assumptions of steady flow, no 
energy transfer to or from the fluid, no body forces, no 
shear stress (inviscid), and irrotational flow the velocity 
potential, <t> , is utilized [Ref. 3:p. 48]. 4>, a scalar 
function of spatial coordinates, x and y, is defined such 
that 

V = V«t» (3.1) 
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and 



04> 


c)4> 


c)x 


V — 

c)y 



The importance of the velocity potential is that only one 
equation is needed to describe the irrotational flow. 
Velocity components are obtained using equation (3.2). 



B . LAPLACE EQUATION 

For steady, incompressible flows the continuity equation 
(2.9) is 



c)u dv 
— + — =0 
c>x dy 

Rewriting (2.9) in terms of the velocity potential the 
equation becomes 

£) 2 4 > 3 2 4 > 

+ = 0 (3.3) 

clx 2 c)y 2 

This form of the Laplace equation [Ref. 2:p. 81] is the 
governing equation for steady, irrotational flow of an 
incompressible fluid. 

The importance of equation (3.3) is that it is linear 
allowing for the principle of superposition. For example if 
<j>i/ 4>a/ 4> 3 . . . are solutions of (3.3), then the sum = 4>i + 

+ 4> 3 +...is also a solution of (3.3). Superposition of 
irrotational, incompressible flows allow for the construction 



13 



of complex flows that are also irrotational and 
incompressible . 

C. SIMULATION (CONFORMAL MAPPING) 

The inviscid flow about an airfoil can be obtained most 
conveniently by means of a transformation, which starts with 
a known flow about a simple contour, a circle, distorts the 
contour into the desired shape, and simultaneously adapts the 
flow to that shape. The transformation is accomplished using 
a sequence of three conformal mappings [Ref. l:p. 324]. 

The first mapping, necessary only when the airfoil 
trailing edge has non-zero thickness, is accomplished using a 
logarithmic mapping function. The airfoil is perturbed 
slightly to make the upper and lower surface, trailing-edge 
points coincide. 

The second mapping analytically removes the trailing-edge 
corner using the Karman-Tref f tz mapping. 

The third and final mapping transforms a quasi-circular 
shape into a perfect circle using an iterated sequence of 
Fast-Fourier Transform applications. 

During the transformation of streamlines about a circle 
to those about an airfoil, the preferable approach insuring a 
transformed flow free from vorticity uses the complex 
variable z = x +iy [Ref. 5:p. 285], The transformation of z 
to another plane is Z= f(z) = £, + ii) . The potential 
function, Q = <J> + i+, is irrotational and incompressible in 
both planes. 
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The streamlines, <l> , and equipotential lines, 4>, of a flow 
in the z-plane will transform into another orthogonal network 
of lines in the C “Plane. Different magnification ratios and 
different amounts of rotation at different points in the 
field will, however, change the appearance of the flow 
pattern from about a circle to about the airfoil. 

The general transformation function whose derivative 
dC/dz satisfies the requirement dC/dz approaches unity for 
large values of z is 

C= Z + Col nZ + Ci/z + C 2 /z 2 +... . (3.4) 

The requirement is necessary as streamlines are not distorted 
a great distance from the body where the body's shape has no 
influence on the flow. 

The coefficients may be real, imaginary or complex. A 
finite number of the coefficients are determined from the 
specified normal velocity components equally spaced around 
the unit circle, and from the Kutta condition which ensures 
stagnation at the trailing edge. 

While in the first iteration the normal velocities are 
zero, and the solution for flow over a circle is used, the 
subsequent normal velocity boundary conditions are determined 
from the previous viscous-flow calculations using the 
equation 
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(3.5) 



d 

V„ = — (u e 6~) 
ds 

where u e is the velocity at the edge of the boundary layer 
and 6* is the displacement thickness. Once the coefficients 
are found, the real and imaginary parts of equation (3.4) are 
equated yielding 



i, = 4(x,y) and t) = t)(x,y) . 

As x 2 + y 2 = r 2 the two equations of E, and t) are transformed 
to 



Then x 2 



and y 
x 2 + 



x = x( E, , r 2 ) and y = y(r),r 2 ). 

2 are added to yield 
y 2 = r 2 = Qx ( E, , r 2 )^j 2 + Qy(r),r 2 ) 




(3.6) 



After dividing both sides by r 2 



1 




1 




r 2 


xU,r 2 ) 


2 + — 
r 2 


y(T) ,r 2 ) 



Then each circle of radius, r, in the z-plane is transformed 
to the proper shape in the C -plane to describe inviscid flow. 
1. Transformation of Velocities [Ref. 5:p. 291] 

In the z-plane as Q(z) = 4>(x,y) + il(x,y) the 

velocities are defined by 
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dQ(z) 



dz 



= V* - iv. 



where , 



V = V x + iv y . 



In the £ -plane the velocities are also found 



( 3 . 8 ) 



dQ 

dt 



= V*. - iV^ 



( 3 . 9 ) 



where , 



V C = v 4 + iV v 



The velocities in the two planes are equated by 



v 4 - iv n = - * - = 



dQ dz Vx - iV y 
dz dC dC/dz 



( 3 . 10 ) 



The pressure in the transformed stream is related to the 
stream velocity through Kelvin's equation 



1 




1 




- p 


v c 


2 + P = Constant = — p 


V C 


2 


2 



+ Pa. 



( 3 . 11 ) 
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IV. VISCOUS FLOW 



A. DERIVATION OF BOUNDARY LAYER EQUATIONS 

The previous analyses provide a valid solution to the 
flow outside the boundary layer. Within the boundary layer 



however. 


the 


effects of 


viscosity cannot 


be 


neglected . 


In 


laminar 


flow 


governing 


equations can 


be 


obtained 


by 



simplifying the conservation equations. In turbulent flow, 
however, the number of variables outnumbers the equations. 
Great dependence is then placed on dimensional reasoning and 
on hypotheses suggested by experimental results. 

The most important deduction from Prandtl's thin boundary 
layer theory is that static pressure can be considered 
constant across the boundary layer [Ref. 3:p. 299]. 

Op 

— = 0 (4.1) 

Oy 

As the boundary layer thickness, 6, is small, d6/dx is also 
small. Streamlines are then only slightly curved and the 
radii of curvature, R, are large. With a large R the 
equilibrium condition 



Op u 2 

Oy P R 
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illustrates that Dp/ Dy will be very small and can be 
neglected. Experimental results confirm that Dp/ Dy may be 
neglected even over surfaces of small radii of curvature. 

Also, in a thin boundary layer with a slowly changing 
thickness Dv/Dx is much smaller than Du/Dy. The significance 
is then that the normal shear stress may be neglected when 
compared with the viscous shearing stress. Equation (2.18) 
then becomes 



Du 

Tjty — P . (4.2) 

Dy 



With this simplification the approximate equation for 
flow in the two-dimensional boundary layer can be found 
directly. Newton's second law of motion applied to a fluid 
element of mass may be written 



Du 


Du 


3u 


C>P C>T V * 




pdxdy ( — 


+ u 


+ V—) 


= (- — + 


• ) dxdy (4.3) 


Dt 


dx 


dy 


Sx c)y 




as illustrated in 


Figure 


4.1. 


Substituting 


equation (4.2) 



and dividing both sides by dxdy yields 



Du Du Du Dp D Du 

p( — + u — + v — ) = ( - — + — (y — ). (4.4) 

at Dx Dy Dx Dy Dy 



In terms of kinematic viscosity equation (4.4) becomes 



Ou 


du 


Du 


i Dp 


D 2 u 


+ 


u + 


V — 


— 


+ V — 


dt 


Sx 


Dy 


p Dx 


Dy 2 
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Pdy — 



o y3t 

(o yJt + dy)dx 

y 



■> 



— > 

dy 

dx 

< 

a yit dx 



A 



y 



> x 



P 

< (P + — dx)dy 

x 



Figure 4.1 Forces Acting on an Element in the Boundary Layer 
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This equation is the boundary layer equation of motion and is 
identical to the equation found using an order-of-magnitude 
analysis [Ref. 3:p. 443]. Equation (4.5) is also nearly 
identical to the Navier-Stokes equation (2.23) with the 
exception that the term vc) 2 u/Bx 2 is deleted. The order-of- 
magnitude analysis suggests that this term, vS 2 u/c)x 2 , may be 
neglected compared to vd 2 u/ dy 2 . Combined with the 
continuity equation c)u/dx + Ov/Oy = 0 (2.9), equations (4.5) 
and (4.1) are known as Prandtl ' s boundary layer equations 
[Ref. 2 : p . P-9]. 

For an incompressible flow, there are three variables, u, 
v and p, but only two equations, (4.5) and (2.9). The 
equations may be solved though, by first determining p as a 
function of x using inviscid methods, setting 0p/c)y = 0 in 
the boundary layer, and then solving (4.5) and (2.9) for the 
velocity distributions. 

B . TURBULENT FLOW 

Turbulent flow as differentiated from laminar flow is 
characterized by fluctuating instantaneous properties which 
greatly increase the complexity of the problem. A very 
useful simplification to the turbulent problem is then the 
use of time-averaged values, denoted by a bar over the value. 
Instantaneous values are indicated by the prime [Ref. 4:p. 
23 ] . 
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u = u + u' 
v = V + V' 
p = p + p' 

The continuity equation containing total values becomes 



0 _ 0 _ 

(U + U') + (V + V') - 0. 

3x c)y 



Simplifying the equation becomes 



3 _ 3 3 _ 3 

— (U) + — ( u ’ ) + — (v) + — (v’) = 0 

3x 3x c)y Qy 



with 




3 _ 3 

— (u ) and — ( u ' ) 

3x 3x 



3 

— (u’ ) = 0 

3x 



3 _ 3 _ 

— (v) = — (V) 

3x 3x 



3 

and — ( v 1 ) 

3x 



3 

= —(V ) = 0. 

3x 



The time-averaged continuity equation for turbulent flow is 
now 



3 _ 3 _ 

— (u) + — (v) = 0 (4.6) 

3x 3y 
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Applying total values, the steady version of the Navier 



Stokes equation (2.23) becomes 



0(u + u' ) _ c)(u + u' ) 

(u + u') + (v + v')- 



Ox 



Oy 



i 0(p + p' ) 

p Ox 



+ V 



0 2 (u + u' ) c) 2 (u + u' ) 

+ 



>X‘ 



0y : 



(4.7) 



After simplifying, the time-averaged Navier-Stokes equation 
for turbulent flow [Ref. l:p. C-10] becomes 



_0 _ 

u— (u) 

Ox 



+ V 



0(u) 

Oy 



l Op 0 2 u 0 2 u 

- — + v ( — — + — ) 

p Ox Ox 2 Oy 2 



0 — 0 

- — (u' 2 ) - — (u' v' ) . (4.8) 

Ox Oy 

The new terms, 0/Ox(u' 2 ) and 0/0y(u'v'), which correspond to 
normal and shear stress, are called Reynolds 
stresses. The similar y-component terms are 0/0y( v ' 2 ) and 
0(v'u' )/0x. 

C . TURBULENCE MODELS 

The time-averaged Navier-Stokes equation is nearly 
identical to the original equation except that the 
instantaneous values are replaced by the mean or time- 



23 



averaged values and two additional terms involving 
fluctuating velocities, u' and v', appear. An interpretation 
of these two terms compares them to the previously existing 
terms 0 2 u/c)x2 and c) 2 u/c)y 2 . The right hand side of equation 
(4.8) less the pressure term, and after multiplying by 
density, becomes 

£> 2 u c) 2 u 3 a — 

y + y - p — u ' 2 - p — u ' v ' 

c)x 2 Qy 2 dx c)y 



or 

c) c)u c) c)u 

— (y— - pu' 2 ) + — (y — - p u ' v 1 ) . 

ox ox c)y c)y 

As each term has the dimensions of stress, and y(c)u/c)y) is 
part of the laminar shear stress x v>t , it appears that the 
term -pu'v' represents a turbulent addition to shear stress 
[Ref. 2:p. T-2]. Now, this shear stress is really a vertical 
mixing of horizontally, travelling fluid particles. A model 
of this mixing then calculates the rate of momentum transfer 
involved. 

1. Prandtl's Mixing-Length Model 

To predict the turbulent stresses Prandtl assumed 
that turbulent fluctuations are primarily the result of the 
mean velocity differences between two layers in the flow. 
Therefore, u' is proportional to c)u/c)y with the 
proportionality factor having the unit of length. 
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(4.9) 



(U’ 2 ) * 




Also, assuming that v' is of the same order of magnitude as 
u' at a particular point, 



c)u 

(v' 2 )* = b — . (4.10) 

dy 



Substituting u' and v' into the turbulent shear stress, x T , 
is 



du c)u 

- pu'v' = - pab( — )( — ). (4.11) 

c)y c)y 

As a and b are both unknown constants of length, they both 
may be replaced by the "mixing length", 1, the hypothetical 
distance between the two layers involved. The turbulent 
shear stress is now 



It 



P i 2 



c)u Qu 

dy c)y 



(4.12) 



which insures the correct sign. 

2. Eddy Viscosity - Cebeci-Smith 

The turbulent addition to shear stress may also be 
modeled in terms of "eddy viscosity". As laminar shear 
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stress = Tyx = p(du/c)y) = pv(c)u/ dy, turbulent shear stress 
may be equated 



c)u 

- pu'v' = pv t — , (4.13) 

c>y 

where v t , the eddy viscosity is empirically determined. 

The method used here is that of Cebeci and Smith as 
modified by Cebeci, Clark, Chang, Halsey and Lee, [Ref. l:p. 
327 ] . 

Eddy viscosity by Cebeci is represented by 



{0 . 4y [ 1 - exp (- -) ]} 2 

A 



T t r (0 < y < y c ) (4.14) 



Vt = 



a 



( u e - u)dy 



Y tr Y (Yc < y < 6) (4.15) 



where 



A 



26 v 




and 



1 

Y = . 

1 + 5 . 5 ( y/6) 6 



The distance from the body y c which is less than the boundary 
layer thickness, 6, is the distance where the two equations 
(4.14 and (4.15 give the same resultant v t . 
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The intermittency factor, if tr , which indicates the 
local fraction of turbulent flow to total flow, is given by 



= 1 


- exp 


- G ( X - X t r- ) 


x dx 






■ 


X t r U a 



The location of the start of transition is x tr , and the 
empirical factor G is 



1 u© 3 

G = R- 1 - 3 * (4.17) 

1200 v 2 



where R xt r is the transition Reynolds number 

Rxtr = (U e X/v )tr. (4.18) 

In equation (4.17) the empirical constant GY tr - = 
1200, was used by Chen and Thyson [Ref. 4:p. 327]. Values 
lower than 1200 may give better results at low Reynolds 
numbers as will be discussed in Section VII. 

The expression for a is 



0.0168 
a = 

p2 . 5 



0.0168 



1 



— i 2.5 



c)u/c)x 



B 

c)u/c)y 



(4.19) 



The non-dimensional factor F represents the ratio of 
the product of the turbulent energy by normal stresses to 



27 



that by shear stress evaluated at the location where shear 



stress is maximum. 



F = 1 



(u ' 2 - v ' 2 )c)u/c)x 



- u'v' c)u/c)y 



(- u' v' ) ; 



(4.20) 



According to the data of Nakayama [Ref. l:p. 327], ft can be 
represented by 



ft = 



U 1 2 — v 1 2 



(- U'v' ), 



Rt ^ 1 



1 + 2R t ( 2 - R t ) 



1 + Rr 



Rt > 1 



where R T = x w /(- u'v')max and x w is the wall (body) shear 
stress . 



D. TRANSITION 

The location of the onset of laminar- to-turbulent 
transition when not found experimentally is determined 
empirically. The method used by Cebeci [Ref. 4:p. 333] is 
the criterion proposed by Michel. 

At the point of transition the Reynolds number based on 
momentum thickness, 0, is related to the Reynolds number 
based on the coordinate position, x. 
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22400 

JWr- = 1.174 (1 + )R. xer 046 (4.21) 



where 



Rex = U e (x/v ) and R otr - i U e (0/v). 



In the Cebeci Code the transition may be determined in 
the following ways: 

1) The points of transition are calculated using Michel's 
criterion . 

2) If laminar separation occurs forward of the criterion 
points, Michel's criterion is disregarded and 
transition is redefined at the separation point. 

3) The transition locations may be specified by the user 
provided stall is not computed. 
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V. VISCOUS METHODS 



Momentum transfer in fluids is accomplished by 

hydrostatic pressure and viscous stresses. When viscous 
stresses are negligible, fluid behavior can be predicted by 
inviscid flow methods as stated in Section III. 

Viscous stresses caused by a variation in velocity in a 
direction normal to the flow are called shear. The most 
common shear is that found in the boundary layer between a 
displayed stream and the solid surface. With the "no slip" 
condition fluid velocity is zero on the surface, but the 
velocity gradient is not so constrained. From the body along 
its normal direction the fluid velocity asymptotically 
approaches that of the free stream. 

As mentioned in Section III, Prandtl hypothesized the 
division of the flowfield into the two regions, the boundary 
layer where viscous effects cannot be neglected and the 
region outside the boundary layer which may be considered 
inviscid . 

This hypothesis allows for the use of the parabolic 
boundary layer equations of section III instead of the 
elliptic Navier-Stokes equations. Depending on the boundary 
conditions, solutions fall into three methods [Ref. 6:p. 13]: 

1) The direct boundary layer method. This method uses the 
"no slip" condition, where normal and tangential 

velocities are zero at the surface, and a 

pre-determined velocity is specified at the boundary 
layer edge. 
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2) The inverse boundary layer method. Boundary conditions 
are replaced by wall shear or displacement thicknesses. 

3) The interactive boundary layer method. The edge 
boundary condition drives a combination of displacement 
thickness and external velocity. 

Methods one and three will be discussed as they apply to 
the Cebeci Code. 

A. DIRECT BOUNDARY LAYER METHOD [Ref. 6:p. 13] 

This method for solving the boundary layer equations is 
used only near the leading edge where the viscous effects are 
small. Initial conditions are generated at the stagnation 
point and the equations are integrated around the leading 
edge. The numerical solution utilizes a finite difference 
method where the continuity and momentum equations are 
redefined as a system of linear algebraic equations. 

The method begins by describing steady, incompressible, 
2-D flows in a curvilinear coordinate system where x is 
directed along the airfoil surface and y is perpendicular to 
x. The boundary layer equations with the turbulent Reynolds 
stress are 



du dv 

— + — =0 (5.1) 

dx dy 



l dp 


l d 


du 


du 


_du 


— + 


— 


[u— - 


pu'v’ ] = u + 


V 


P dx 


P dy 


dy 


dx 


dy 
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(5.3) 



dp 

— = o. 
dy 



where the order of magnitude of any turbulent stress is 
assumed to be that of its laminar stress. The boundary 
conditions are: 



aty=0;u=0,v=0 (5.4) 



at y = oo ; u = Ue ( x ) (5.5) 



The eddy viscous stress previously defined is reprinted as 



du 

-pu'v' = pv t — . (4.13) 

dy 



Also, the pressure gradient term may be written 



1 dp du e 

— — = u e - — (5.6) 

p dx dx 



Therefore, the momentum equation (5.2) may be rewritten as 



_du 


+ du 


du e 


d 


du 




u 


V 


U e + 





(b — ) 


(5.7) 


c)x 


dy 


c)x 


dx 


dy 



where b = 1 + v t /v and the boundary conditions are 

at y = 0; u(x,0) = 0, v(x,0) = 0 (5.8) 

at y = y e ; u(x,y e ) = u e (x) . 
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1. Falkner-Skan Transformation [Ref. 6:p. 14] 

To solve the new boundary layer equations the 
Falkner-Skan transformations are used, which reduce the 
number of variables, and scale the normal coordinate y and 
the stream function 4* with reference to the external 
velocity. 



(5.9) 

<Wx,y) = (u. x)* ♦ f(x,n) (5.10) 

The continuity equation is automatically satisfied 
using the stream function (u = O'I'/Oy and v = -0^/Ox) . 
Therefore, only the momentum equation needs to be solved, 
which after transformation becomes 

m + 1 Of ’ Of 

(bf")’ + ff" + m [ 1 -(f) 2 ] = x(f* — - f" — ) ( 5 . 1 1 ) 

2 Ox Ox 

where m = (x/u e ) (du e /0x) , a dimensionless pressure-gradient, 
and f' = Of /On • 

This equation (5.11) is a third order partial 
differential equation, and the solutions are "non-similar" as 
they are functions of both x and n- If the solutions were 

only a function of n > then the right hand side of the 
equation would equal zero, and the flow would be "similar" 
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[Ref. l:p. v-10]. To solve equation (5.11), a numerical 
solution is needed such as the "box" method. 

2 . The Box Method 

The box method, developed by Keller in 1970 [Ref. 
l:p. 331], is a widely used methods for solving non-linear 
differential equations. The steps of this method include the 
conversion of the Falkner-Skan , transformed, momentum 
equation into a system of first-order partial differential 
equations. This non-linear system, after conversion from a 
continuous function to discrete, is linearized by Newton's 
method. The block elimination method is then used to solve 
the linearized difference equations of the boundary layer 
problem . 

The third order momentum equation (5.11) is converted 
into a first order system with the addition of the dependent 
variables U and V [Ref. 6:p. 14]. 



U = f’ 



(5.12) 



V = U' = f" 



(5.13) 



m + 1 Du Of 

(bV) ' + fV + m ( 1 - U 2 ) = x(U — - V — ) 

2 Ox Ox 



(5.14) 



The boundary conditions are 

at i) = 0; U(x,0) = 0, f(x,0) = 0 
at r) = I) U ( X , T] e ) = 1 
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The solution domain, 0<x<x x and 0<r|<r)j, of the 
continuous functions f, U and V is then covered by a 
rectangular grid to facilitate the problem solving with a set 
of discrete values. This grid is shown in Figure 5.1. 

The subsequent notation [ is used to represent the 
quantities of f, U or V at the point (x*, t^j). 

All quantities can then be approximated by network 
coordinate values. Using the grid system, the solution of 
the parabolic layer equation at a certain streamline position 
depends solely on the solution of upstream positions, while 
no downstream influence needs to be considered. As 
calculations proceed from the stagnation point in the 
downstream direction, the overall solution can be obtained 
incrementally. Hence, one step of the solution procedure 
sets up the governing equations for a column of grid boxes in 
the sub-domain 



Xi-i<x<Xi and 0 <ti<tij 

and solves for the values of the downstream grid position. 
The x-grid position currently solved for is then assigned the 
superscript "i" while "i-1" represents the previous position 
of known properties. Using coordinates of box midpoints and 
centered-dif f erence derivatives, the equations are actually 
satisfied midway between the grids. 

Equations (5.12) through (5.14) in terms of finite 
differences [Ref. 6:p. 15] are now written 



35 




LEGEND 

O known 

□ unknown 

A center for momentum 
equation 

V center for equations 
containmq 
^-derivatives only 



Figure 5 . 1 



Rectangular Grid for Finite Difference 
Approximation 
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(5.15) 



f 1 - f 1 
1 j-l 

“ *5 



7 (u 3 + U J-1> 



U j - U J-l 



7 (v d + v j-!> 



< bv >r 4 - (bV, j-l . ♦ 1 

. + ^ 






i-l 



(fv); i + 

J~2 



m 



i - ! E- (v ^>3 = 



X 



1-2 



u 1 - u 1 ' 1 

tt i- 2 U j-2 U j-2 

U j-| TT< 



(5.16) 






(5.17) 



where the ordinary differential equations, (5.15) and (5.16), 
are centered about (x±, rjj-l) and the partial differential 
equation, (5.17), is centered about th-*). 

The boundary conditions are 

U* = 0, fij = 0 and Uj = 1 

Equations (5.15) and (5.16) are the centered difference 
derivatives . 

3. Newton's Method [Ref. 6:p. 15] 

This set of finite difference equations is nonlinear 
with combinations of unknowns. Newton's method is therefore 
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needed to solve the system. The variables are linearized 
using an iterative procedure with preceding values. 



f*' 1 = fj _1 + 6f j ' 1 and f* ,K 
where, 6fJ ,K << fJ ,K_1 
U*' 1 = Uj _1 + &U*' 1 and U*' K 
where, 6U^ ,K << Uj ,K_1 

V j 1 = V j~ 1+ 6v j 1 and V j' K 
where, 6vV K << Vj ,K_1 



fj ,K-1 + 6fj ,K for K>2 



U*' K_1 + 6j' K for K>2 



V*' K_1 + 6V* ,K for K>2 



K is the iteration counter. After substituting these values 
into equations (5.15) through (5.17) and neglecting the 

i , K i , K i.K 

quadratic terms of 6f , 6U and 6V , the system of 

j j 3 

unknowns is then linear. 



c^i.K c C i-K 

6£ j ■ 6f j-i 

+ h., u *;?” 1 

3 3 s 

6uV K - 

+ h 

3 3 2 



h . 



— ( 6U^ ’ K + 6U*'!f) 
3 3~1 



= f 



i.K-1 

j-1 



- f 



i.K-1 



_^2( 6V^ • K 



+ sv^:^) 



"J:?- 1 



i.K-1 

3 



(5.18) 



(5.19) 
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( S a ) j ’ K 6V j ' K + ( S 2 ) j * K 6V j ' K + (S 3 )J ,K 6fj ,K + 



i . K c -i . K 



i . K r „i . K 



i.K„wi.K 



(S*) j 6f y_“ + (S s ) + ( S 6 ) j 6U j _ x 



(r 2 )j‘ K 



where 



b 1 ^’ 1 x 

/c , i.K_ c j A x i-| /JC i.K-l 

S 1 j FT - + 






i-l m 1+1 i.K- 

D-r + f j-i 



(S 2 ) j- K = 



b 1 ^" 1 X 

3 + X i*i {f i.K-l 

hj + (f j-J 



+ ^T 1 f j 



(sa) j ,K= tk~ ^i-T 1 + v j : !> + v j’ K_1 

(s ^j' K = 7 F 2 tVjlj ' 1 + v j-^ + v j l ?" 1 



( S5 ) 



i.K 



- (fjli + 

K i 11 



(s 6 ) J ,K = 



/ X i — 2 , i UT i . K-l 

( “T — + m )U j-l 



/ .i.K 
( r 2 ) j = 



/ V>T r \ 1 " 1 

(kV) . ~ m + 1 1 K— 1 

{( 1 j- LJ_ + !L_^_L(fv)^ 



X 



i - 2 . _i . ,,.i . K-l 



x i - 2 ,.,i . K-l ^i.K-1 



- <-^ + > 2 + q: 



(5.20) 



1 



• K-l 
-1 
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+ f ^:^ _1 

3 2 5lJ 

(bv)* 

- {( i- 






3-2 



i-2 



_ 1 

2 



- (bvjj., 



m 1_1 + 1 






3 2 



+ (• 



i-| . _i-l, ,„i-l 



T - 

i 



+ m 



) (u: 



3-2 






+ 2m 1 -*} 



The boundary conditions are 

6f o = 0, 6U 0 = 0, 6Uj = 0 
4. Block Elimination Method [Ref. 6:P. 17] 

The system of equations are iteratively solved until 
i.K i.K i.K 

6f , 6U and 6V become small enough to be neglected. 

3 3 3 

The solution method is that by Keller and is called the Block 
Elimination Method. In this method block-tridiagonal 

matrices are composed of blocks. Only those blocks on the 
main diagonal and on the two ad3acent diagonals are non-zero. 
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[Aj' K ] 


[c i 


k i 














<&i’ K y 


< r i' 


K > 


[B 2* K] 


1 A 2 


k i 


t C 2 


K i 












(r ': 


K > 




ib i; 


k i 




K i 




K] 






{6 j ’ K } 


■ <r ^: 


K }(5 










K i 
■1 J 




K J 
-1 1 


[C J- 


K ] 
-1 J 


/ c i • K > 


< r j: 


K 

■l t 












1b j' 


K J 


t A j' 


K ] 


{6j‘ K > 




K > 



where the blocks are 3x3 matrices 





T 


-hj/2 


0 


[ A f ’ K ] = 


(S 3 )j - K 


(S s )j- K 


( S i ) j ’ K 




0 


-1 






“1 


-hj/2 


0 — 


[b* , k ] = 


(S.)j' K 


(S.) j ,K 


(Sa) j ,K 




0 


0 


0 




u 


0 


0 — 


[Cj ,K ] = 


0 


0 


0 




0 


1 


" h j+l /2 




T 


0 


0 — 


[A^- K ] = 


0 


1 


0 




0 


-1 


-V 2 _ 



2< j < J-l 



2 < j < J 
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[aJ- K ] = 



-V 2 



(S3) 



i.K 



(Sb) 



i . K 



(Sa) 



1 . K 



The right-hand sides of the equations are 



( r* , V K = + h.U 1 *?- 1 

3 3-1 3 3 3-2 



2<h< J 



(r 2 )^' K = {As defined in equation (5.20)} 2<j<J 



, .i.K 
(r 3)j 



i.K-i _ yi-K-1 h V 1.K-1 
3 3-1 D+l 3+2 



2<h< J 



( r ) ^ ’ K = 0, (r 2 )i ,K = 0, (r 3 ) j‘ K = 0 



The block elimination method solves the linear 
equations with two steps. The forward step eliminates the 
lower diagonal of the tridiagonal matrix. The reverse step 
solves the remaining system from bottom to top. 

B. INTERACTIVE BOUNDARY LAYER METHOD [Ref. 6:p. 18] 

As the direct boundary layer method is, as previously 
stated, restricted to regions of small viscous effects, and 
integration of the boundary layer equations fails at points 
of zero skin friction, a method is needed to integrate the 
boundary layer through the point of emerging reversed flow. 
This method must also account for strong interaction effects 
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due to separation and rapid acceleration of the flow 
downstream of the trailing edge. 

The interactive method fulfills these requirements by- 
treating the external velocity and displacement thickness as 
unknown quantities. Reflecting the elliptic nature of the 
outer flows, an additional unknown is introduced, but the 
solution can be obtained using either the eigenvalue or 
Mechul function methods. 

The Mechul method is preferred as the eigenvalue method 
involves nonlinear problems. In this method the edge 
boundary condition of the direct method is supplemented with 
the interactive boundary condition. The unknown external 
velocity is related to its displaced and perturbed 
conditions. The unknown functions u(x,y) , v(x,y) and u e (x,y) 
are represented in this system of boundary layer equations 



Ou 


3v 







-f- 


= 0 


dx 


c)y 




c)u C>U 


c)u e 


c> c)u 


U + V = 


u e 


+ v — (b — ) 


c)x c)y 


c)x 


c)y c)y 



c)u e 



= 0 

3y 



( 5 . 22 ) 



( 5 . 23 ) 

( 5 . 24 ) 



where pressure in the y-momentum equation is expressed in 
terms of the external velocity. 

The Mechul function approach assumes that the external 
velocity, u e , is a function of two arguments, x and y, 
allowing for an easy setup of finite difference equations, 
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and avoiding nonlinear eigenvalue techniques. The velocity 
components u and v are required to satisfy the no-slip 
condition at the surface, and u must merge smoothly into the 
edge velocity. 



at y = 0 : u(x,0) - 0, v(x,0) = 0 



at y = y e : u(x, Y e ) = u e (x, Y e ) , u e (x,y e ) = 



= U eT 




d d4 

— ( u e 6* ) 

cH x-4 



U e i(x) is the inviscid edge velocity and the last term, 
the Hilbert integral, approximates the viscosity induced, 
perturbation velocity. 

Interactive methods are useful in both attached and 
separated regions, while direct methods fail at the onset of 
reversed flow, and inverse methods converge poorly. Only at 
the stagnation point singularity are interactive methods 
prohibited. The transformation of the partial differential 
equations into a linear system of algebraic equations is very 
similar to that of the direct method. The normal coordinate 
y, streamf unction <|> , and the external velocity u e are scaled 
with reference to a constant velocity u 0 , and the local 
streamwise coordinate x. 

( 5 . 25 ) 
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<Hx,y) = (u 0 vx)* f ( x , i) ) 



(5.26) 



u e (x f y) 

W ( x , t) ) = (5.27) 

Uo 

Uo is the vector mean velocity. U e cannot be the reference 
velocity as for Falkner-Skan variables because in this case 
the external velocity is unknown. The first order, semi- 
transformed coordinate system with additional variables U and 
V is 



U = f' 




(5. 


.28) 


v = u' = f" 




(5, 


.29) 


1 Ow Ou 


Of 






+ - fv + xw — = X(U — - 


V—) 


(5. 


.30) 


2 Ox Ox 


Ox 


o 

II 

£ 




(5. 


.31) 



and the boundary conditions are 



at i)= 0: U ( x , 0 ) = 0, f(x,0) = 0 

at T) = t) e : U ( x , T) e ) = W(X,T} e ) 



W(x,n«) 



Uel(X) 



Uo 




{( — ) *[W(^,TJe)n« 

Uo 



- fu,n-)]>- 
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The conversion of the flowfield to discrete values is 



very similar to that of the direct method with an orthogonal 
grid, central differences, and two-point averages. In the 
interactive method though, the solution proceeds in the 
downstream direction only. As only downstream disturbances 
are accounted for, backflow causes numerical instabilities. 
A stable integration can be accomplished though, with the 
assumption that backflow velocities are comparatively small. 
The FLARE approximation (Flugge-Lotz and Reyhner) , [Ref. 6:p. 
19] sets the streamwise convection term uc)u/c)x equal to zero 
in regions of backflow. 



The finite difference equations of the interactive 
boundary layer with the "on-off switch" are now 




0 



1 



if U 1 ! >0 
3 ~ 2 

if uj j. <0 

J “ 5 




(5.32) 




(5.33) 



<bV)*“* -(bV)J:| 
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= X 



i-i 



E 



FLR j _ | ( U 1 ” ^ 



U 



- U 
IT 



i-1 



j-i 



- (V 



i-l f 



- f 1 ’’ 



j-5 



W j " ”^-1 



= 0 



(5.34) and (5.35) 



The boundary conditions are also expressed in terms of 
grid or nodal values. A panel method type approximation 
leads to 



= 0 f = 0 









W J' W J = g i 


+ - 




where g± 


and 


c ± ± 


represent a 


parameter 


and 


element 


of 


the 


interaction 


matrix, due 


to 



approximation to the Hilbert integral. To keep the number of 
generated terms to a minimum, ordinary differential equations 
like the y-momentum equation are centered about the 
downstream face, and partial differential equations like the 
x-momentum equation are centered about the middle of the box. 

The unknowns occur in vectors of four components 



{fj, Uj , V*, wj}. 

The J quadruplets of unknowns match with 4J equations, 

including 2(J-1) auxiliary relations and (j-1) x-momentum and 

(j-1) y-momentum equations. Each equation corresponds to one 
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of the (j-1) grid rectangles and four boundary conditions. 
The system is linearized around the values of the preceding 
iteration (counter K-l) and a system with Newton iterates 



c .pi • K ATT i.K c, r i.K At7 i- K 

ofj , oUj , oVj , oWj , emerges. 



6f 



i.K cpi.K h j. KTI i.K - i.K. _ ei.K-1 -i.K-1 

j ■ 6f j-i — r- (6U j + 6a j-i } - f j-i ■ f j 



+ h 

J J 2 



(5.36) 



6U ^ ' K - SuV* 



-ji(6Vj ,K + 6vJ^) 



uj.K-i _ yi.K-1 



+ h . v^;?- 1 

3 J 2 



(5.37) 



(Sx)V K 6vV K + (S a )V K 6vV? + (S 3 ) 1,K 6f 1,K + 
3 3 3 3-1 3 3 



/ P \ 1 ' Kepi » K j e i i • Kerri . K , / e \ i * Kerri * K , 

(S 4 )j 6 f j-i + ( Ss )j 6U j + ( s «)j 6 U j-i + 



/ C \ ^ ' K C tri »K r e \ 1 * Kerri . K er ri «K / \ 1 . K 

(S 7 )j OWj + ( S a ) j 6 W j-l 6 W j-l ( r z)j 



(5.38) 



eeri.K e r , i • K „i-K-l ..i.K-1 

6w j • 6w j-l = w j-l - w j 



(5.39) 



where 



b 1 ^ -1 x 

fs ) i ' K - j + i_ 2 

(Sl) j TT— + 7T~ {f j-l 



i-l 1 i.K-1 

f j-i } + T f j 
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i . K-l 



(S 2 ) j ,K = 



x . 



3 + ” 1 "2 (f 1 - 

“TTJ— + 1Y— (f j- 



i.K-1 



fi-i 1 -i.K-1 

f d-r ^ f j-i 



(S 3 ) j ,K = 



X 



i-4 (yi-K- 1 + yi- 1 ) + 1 yi-K- 1 

21c— (v :-i + V j-S' + 7 v j 



i-K. . x i-2 ^ „i-l, ^ 1 ,, i.K-1 



■ * TST- 1 V j-I + V M> + T 



V 



j-1 



(Ss)J* K = 



X • 



4i_ uH^flr* , 



(Ss) j‘ K = 



X . 



jll FLR* i 

J-1 J-3 



i.K_ X i-| i.K-1 

(S7) j " He— W j 



i.K_ X i-i W i-K-1 
(Sa) j ' He— W j-1 



i K i ” i— 1 1 iK— 1 

(r*)j = - <( 



+ (^i nw 2 )^:^ -1 - ( u 2 )^:^‘ 1 flr ^_ 1 

J 2 J 2 J 2 



♦ vH _1 fj :?- 1 + vj:; tj :?- 1 

J 3 J 3 J 3 J 3 
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The boundry conditions are 



6uJ* K = 0, 6f* ,K = 0, 6Uj ' K 



6W^ ‘ K = 0 
J 




1 




i . K-l 



The overall procedure is a repetitive, linear approach to 
solve the nonlinear system. The numerical solution is again 
obtained using the block elimination method by Keller, except 
that unlike the direct method the vectors of the unknown 
Newton iterates are four dimensional. 




i . K 




,i . K 




and 
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In matrix-vector form 



[Aj- K ] [cj‘ K ] 








{6j* K > 


<rj K > 


[B^ K ] [A^ K ] 


[c*- K ] 






j c i • K > 
{6 2 > 


<r^- K > 


[bJ’ k ] 


[Aj ,K ] [cj- K ] 






{5 j ' K } = 








[Bj:*] [aJ : *] 
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<■£?> 


— 




[ B j* K ] 


[Aj- K ]_ 




<aj-, 


(rj' K » 


The 4x4 matrix 


blocks are 












T 


-h_./2 


0 




0 




[ Aj ' K ] = 


(s 3 );j 


• K (S 5 )j- K 


(Si)j' 


K 


(S 7 ) j- K 


2< j < J- 1 




0 


-1 


- h j + l/ 2 








0 


0 


0 




0 






“1 


-hj/2 


0 




-1 ~ 




[bJ' k ] = 


(S 4 ) j 


• K (S 6 )j* K 


(S 2 )*‘ 


K 


(S e )J- K 


2< j < J 




0 


0 


0 




0 






_o 


0 


0 




0 






U 


0 


0 




0 




[cj' K ] = 


0 


0 


0 




0 


l<j<J-l 


0 


1 


- h :i + l /2 


0 






0 


0 


0 




1 
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[A^* K ] 



IT 

0 

0 

0 



0 
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-1 

0 



0 

0 

-h 2 /2 



0 

0 

0 

-1 



[A^' K ] 
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i . K 



-h /2 
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i . K 
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i.K r i . K-l -i.K-1 



5- 



- f A 



i.K-1 



+ h.V^_j 

3 3 s 



2 < j < J 



(r 2 )V K = {As defined in equation (5.38)} 2<j<J 



( r 3 ) j ’ K = ^' K_1 - + h j+ i v j;f " 1 



l<j<J-l 



(r«) 



i.K 



wi.K-i 

3 



V^i.K-i 
3 + 1 



1< j<J-l 



(roj ,K = 0, (r 2 )J- K = 0 
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g -c f^"^~^-(l 
i.K 9 i ii z J [l 



C )W^ " 

c ii-J ,w J 



c ii 



and 



(r*)j- K = 0 
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VI. INTERACTION METHODS [Ref. 7:p. 79] 



Interactive methods couple viscous and inviscid flows and 
are intended to compute through regions of flow separation. 
Given their levels of success, these methods have become 
inexpensive alternatives to the Navier-Stokes solvers. 

The simple, classical method of computing the viscous 
flow over an airfoil is: 

1) The velocity distribution is computed for inviscid 
f low . 

2) The inviscid velocity is input to the viscous flow. 

3) The viscous flow is computed by integrating the 
boundary layer equations. 

Now, this method is good at predicting lift and drag, but 
only if the flow remains attached, as information is 
transferred only once from inviscid to viscous regions. For 
more complex flow multiple information transfers are 

required. 

Close coupling is needed to compute flows with separation 
or separation bubbles. A better method than the previously 
outlined classical method for exchanging information between 
viscous and inviscid regions is interaction. The different 
elements of interaction include direct and inverse, inviscid 
and viscous flow solvers. Table 6.1 illustrates the 
different elements. 

The disadvantage of the direct boundary layer method is 
that the equations become singular at the point of 
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separation. The point of separation may be integrated 
through, however, if the external velocity is computed with a 
predetermined displacement thickness. This method is known 
as the inverse boundary layer method. 

Another problem associated with separation is the 
instability of numerical methods which prohibits downstream 
marching in regions of reversed flow. In the situation where 
flow is reversed the FLARE approximation is used, where the 
momentum transport term uDu/c)x is neglected. This 
approximation is not necessarily accurate, but it does allow 
for continued calculations. 

Four interaction models have been developed to calculate 
combined inviscid and viscous flows. All procedures solve 
the Laplace equation for inviscid flow and the boundary layer 
equations for viscous flow. The four models are the direct, 
inverse, semi-inverse and viscous-inviscid interaction 
methods. Each model is subject to different boundary 
conditions . 

The first three models are considered weak interaction 
methods in that they provide only a loose coupling between 
viscous and inviscid regions. The two regions are treated 
alternately. As indicated in Table 6.1, the viscous flow 
solver calculates the flow in the viscous region and produces 
the boundary condition of the inviscid region. The inviscid 
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TABLE 6.1 INTERACTION ELEMENTS 





BOUNDARY CONDITION 


Flow 


Direct 


Inverse 


Inviscid 


* Zero normal 
velocity at 
the surface 


* Prescription of 

velocity distribution 


Viscous 


* No slip 
condition 


* No slip condition 




* Prescription of 
external velocity 


* Prescription of dis- 
placement thickness 
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flow solver calculates the flow in the inviscid region and 
produces the boundary condition of the viscous region. The 
weak interaction methods process either displacement 
thickness or external velocity as an input and the other 
quantity as an output. 

In contrast, the fourth method, the simultaneous 
interaction method, is considered a strong interaction 
method. A strong method calculates displacement thickness 
and external velocity simultaneously. The foundation of the 
four interaction methods are discussed below. 

A. DIRECT INTERACTION METHOD 

The direct interaction model is composed of direct 
inviscid and viscous flow solvers as indicated in Figure 
6.1a. The external velocity distribution is calculated first 
by inviscid computations. The displacement thickness, 6*, is 
then calculated using the external velocity as a boundary 
condition. An updated shape of the displacement body is then 
computed, and all steps are recomputed in order until the 
results converge. As previously stated this method breaks 
down at the point of separation, and is therefore not useable 
in regions of separated flow. However, it is very useful 
where viscous effects are small. The direct method is used 
in the Cebeci Code around the nose and stagnation point of 
the airfoil. 
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a) 



DIRECT 



b) INVERSE 



Airfoil 

geometry 




c) SEMI-INVERSE 



Figure 6.1 Organization of Interaction Methods 

a) Direct, b) Inverse and c) Semi-inverse 
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B. INVERSE INTERACTION METHOD 



This method was developed to circumvent the singularity 
problems near separation. According to Figure 6.1b it uses 
inverse inviscid and viscous flow solvers. Because of the 
inverse method's very slow convergence, though, it is 
suitable only at singular points. 

C. SEMI-INVERSE INTERACTION METHOD 

The semi-inverse interaction method incorporates direct 
inviscid and inverse viscous flow solvers such that 
displacement thickness is input to both solvers as shown in 
Figure 6.1c. External velocity is output from both solvers. 
Convergence is ensured with the use of a relaxation formula 
which redefines the displacement thickness distribution. 

U ev (x) 

6new* ( X ) = 6oldl*(X)[l + U ( ~ 1) (6.1) 

Uex(x) 

where w is a relaxation parameter. 

The numerical weaknesses of the direct and inverse 
methods are improved, but inviscid and viscous regions are 
still loosely coupled. 

D. VISCOUS-INVISCID INTERACTION METHOD 

The viscous-inviscid interaction method ensures a strong 
interaction between the outer, inviscid, and inner, viscous, 
regions. Both the external velocity, u e (x), and displacement 
thickness, 6*, are unknown quantities. Convergence is 
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ensured through the interaction law which uses the blowing 
velocity concept. 

The equations are solved through successive sweeps over 
the airfoil surface as indicated in Figure 6.2. For each 
sweep the external velocity for the boundary layer equation 
is written 

u® ( x ) = u eX (x) + 6u e (x) (6.2) 



where , 



u eI (x) is the inviscid velocity 



and 



u e (x) is the perturbation velocity due to the 
boundary layer displacement. 

The perturbation velocity is modeled by the interaction 
law with the help of blowing velocities. The displacement 
effect of a boundary layer is obtained by ejecting fluid at 
the surface of the airfoil as shown in Figure 6.3. 

With a properly arranged blowing velocity source 
distribution on the airfoil surface, the virtual displacement 
body becomes a streamline. 

In determining the source strengths, the displacement- 
body tangential flow condition is represented by 



v(x,6*) d6* 

u a ( x ) dx 



(6.3) 
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(start) 

l 



INPUT of 

1. Airfoil geometry 
2. Angle of attack 
3. Reynolds number 



V w «. <5* 



INVISCID FLOW SOLVER 

Laolece equation it solved by a 
conformal maoping. 

INPUT; Airfoil geometry and angle of attack 
displacement thickness distribution 
and blowing velocity distribution 
OUTPUT: Velocity distribution u*i on 
displacement body 





r u#l 


DATA INT 

Switch from conformc 
Interpolate external v< 


ERFACE 

il mapping to b.l. grid, 
elocity at b.l. points. 



DATA INTERFACE 

1. RELAXATION of the 

product term * 

2. Comoute BLOWING VELOCITY 
distribution V* t 

3. Switch from b.l. grid to 
conformal mapomg grid. 
Interpolate displacement 
thickness distribution end 
blowing velocity distribution 
at boundary points 



U 0 | 



u.. (5* 



VISCOUS FLOW SOLVER 

Boundary leyer equations are solved by a finite difference method. 



DIRECT B.L METHOD 

(Used close tp lead, edge) 

INPUT: External velocity, 

initial guess of 
velocity profile 

OUTPUT: Velocity profile and 
b.l. characteristics 
like displacement 
thickness <5* end 
skin friction 



INTERACTIVE B.L METHOD 

(Uaed anywhere else} 

Both external velocity and displacement 
thickness are treated es unknowns. 

INPUT: External velocity distribution 

end displacement thickness 
diatribution of current (upstream) 
end previous (downstream) cycle, 
initial guess of velocity profile 

OUTPUT: Velocity profile including the 
'viscous* axtarnel velocity u, 
end b.l. characteristics like diapl. 
thickness <3* end skin friction 



OUTPUT of 

1. Pressure distribution 
2. 8.1. characteristics for both 
upper and lower surface 
3. Velocity profiles if required 



( STOP ) 




Figure 6.2 viscid/Inviscid Interaction Method 
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Blade surface along which 
sources are distributed 



Figure 6.3 The Blowing Velocity Concept 
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Model simplifications are achieved through the use of 
these thin airfoil approximations: 

1) The u-velocity component is considered invariant across 
the boundary layer, as the displacement thickness is 
thin enough to consider the differences negligible. 

2) The blowing velocity, v(x,0), is one half the source 
strength, a(x), on an airfoil represented by a straight 
line . 



a(x) 

= (x,0) 

2 



= v ( x , 6* ) 



'6~ 6v 

— dy 
. o 6y 



d6* diu 

= u e + 6* 

6x dx 



a(x) d 

= — (u e 6* ) 

2 dx 



(6.4) 



where ( d/dx ) ( u Q 6* ) is the blowing velocity. 

The blowing velocity once obtained from the source 
strength is then related to the perturbation velocity, 6u e , 
through the use of the Hilbert integral. 



6U e 



l r~ b 0(4) 

— d 

2 II , xn X~4 



(6.5) 



After substituting equations (6.4) and (6.5) into (6.2) 
the interaction law is obtained. 
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fXb 



( 6 . 6 ) 



u e (x) 



u eI (x) 



1 

+ — 
n; J 



d df, 

— (u e 6*) 

<U x-4 



The numerical implementation of the interaction law 
requires some discrete approximation of the thin airfoil 
integral, equation (6.6). Similar to the panel method, a 
piecewise approximation of the continuous blowing velocity 
d(u e 6*)/dx allows for piecewise analytical integration. 
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VII. AIRFOIL STUDIES 



Cebeci's interactive aircode was applied to four airfoils 
over a wide range of Reynolds numbers and angles of attack. 
The computer program results were then compared to reported 
experimental data. Unless otherwise stated, 20 iterations 
were used for each computer run, and laminar-to- turbulent 
transitions were determined internal to the program. The 
significance of the number of iterations will be discussed 
later in the section. 

A. NACA 66 3 -018 

Computer results of the NACA 66 3 -018 airfoil section were 
compared to the test results of Gault [Ref. 8], which were 
performed in the NASA Ames Research Center 7-by-10-foot wind 
tunnel. The laminated pine model with a 1/8 inch-thick 
mahogany plywood veneer spanned the 7-foot dimension to 
simulate two-dimensional flow. 

Total-and static-pressure surveys, hot-wire-anemometer 
observations, and detailed pressure-distribution and liquid- 
film measurements were made in regions of separated flow. 
The measurements were obtained for a wide range of angles of 
attack and for Reynolds numbers from 1.5 to 10 million. A 
main purpose of these measurements was to identify locations 
of separation, transition and reattachment. 
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Using the Cebeci Code the 66 3 -018 airfoil shown in Figure 
7.1 was initially tested for section lift coefficients. 
Comparisons were made with Abbott and Doenhoff [Ref. 9] and 
the results are presented in Figures 7.2 and 7.3 for Reynolds 
numbers of 3 and 6 million, respectively. 

Upper surface, laminar to turbulent, transition locations 
are shown in Figures 7.4 and 7.5 for increasing angles of 
attack and for Reynolds numbers of 3 and 6 million, 
respectively. Gault's locations were obtained from pressure 
and hot-wire measurements, which provided near identical 
results. The program transition locations were computed to 
be the point of laminar separation. Note that the transition 
locations shift forward as the angle of attack is increased, 
and they approach the leading edge above 6 degrees. Unless 
otherwise stated, all computer runs used a transition 
constant of GY = 1200. 

Midchord upper surface transitions at less than two 
degrees angle of attack and Reynolds numbers of 1.5 to 10 
million are shown in Figures 7.6, 7.7, 7.8 and 7.9. In all 
cases the computed predictions were forward of Gault's 
because of laminar separation predicted by the Code. 

While Gault found leading edge separation bubbles, the 
Cebeci Code did not predict them at any angle of attack for 
Reynolds numbers of 3 and 6 million. 

The relationships between separation and transition are 
illustrated in Figure 7.10 for the results of Gault and the 
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Figure 7.2 Lift Coefficient, NACA 66 3 -018, R = 3 Million 



67 



0.0 0.1 0.2 0.3 



LIFT COEFFICIENT 







Figure 7.3 Lift Coefficient, NACA 66 3 -018, R = 6 Million 
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Figure 7.4 



Upper Surface Laminar to Turbulent 
Transition, NACA 66 3 -018, R = 3 Million 
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LAMINAR TO TURBULENT TRANSITION 




Figure 7.5 Upper Surface Laminar to Turbulent 

Transition, NACA 56 3 _ 013, R = 6 Million 
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Figure 7.6 Upper Surface Transition, Midchord, NACA 
663 -OI 8 , R = 1.5 and 2 Million 
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Figure 7.7 Upper Surface Transition, Midchord, NACA 
66 3 -018 , R = 3 and 4 Million 
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Figure 7.8 Upper Surface Transition, Midchord, NACA 
663 -OI 8 , R = 6 and 8 Million 
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Figure 7 . 9 



Upper Surface Transition, Midchord, NACA 
663 -OI 8 , R = 10 Million 
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Figure 7.10 Upper Surface Transition and separation. 
Midchord, NACA 66 3 -018, R = 2 Million 
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Cebeci code. The experimental results show separation prior 
to transition, whereas the computer results predict 
separation after transition. The importance of this 
difference manifests itself in the difference between the 
measured and computed midchord bubbles, pressure 
distributions and velocity profiles shown in Figures 7.11 to 
7.31. 

Midchord separation bubbles for angles of attack of 0 and 
2 degrees and Reynolds number of 2 million are plotted in 
Figures 7.11 and 7.12. The lines represent contours where 
u / Ue = 0. The Cebeci Code midchord bubbles are much smaller 
than those found by Gault. 

Full chord pressure distributions for angles of attack of 
zero and two degrees, and Reynolds numbers of three and six 
million are shown in Figures 7.13 to 7.16. In each case the 
biggest difference between the experimental results and the 
Cebeci code occurred near the midchord separation bubble 
regions. Figure 7.17 shows a leading edge pressure 
distribution for an angle of attack of six degrees and a 
Reynolds number of three million. 

Midchord velocity profiles are shown in Figures 7.18 
through 7.24 for an angle of attack of zero degrees and a 
Reynolds number of two million, and in Figures 7.25 through 
7.31 for two degrees angle of attack and a Reynolds number of 
two million. These velocity profiles clearly show a big 
difference in bubble sizing. 



76 




Figure 7.11 Midchord Bubble Shape, NACA 66 3 -018, 
AOA =0°, R = 2 Million 




Figure 7.12 Midchord Bubble Shape, NACA 66 3 -018, 
AOA = 2°, R = 2 Million 
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Figure 7.13 Upper Surface Pressure Distribution, NACA 
66 3 -018, AOA = 0°, R = 3 Million 
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Figure 7.14 Upper Surface Pressure Distribution, NACA 
663-OI8, AOA = 0°, R = 6 Million 
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Figure 7.15 Upper Surface Pressure Distribution, NACA 
663-OI8, AOA = 2°, R = 3 Million 
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Figure 7.16 Upper Surface Pressure Distribution, NACA 
663 -OI 8 , AOA = 2°, R = 6 Million 



81 



PRESSURE COEFFICIENT 



UPPER SURFACE PRESSURE 




Figure 7.17 Leading Edge Upper Surface Pressure 

Distribution, NACA 66 3 -018, AOA = 6°, 
R = 3 Million 
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Figure 7.18 Upper Surface Velocity Profile, NACA 66 3 -018, 
X/C = .60, AOA = 0°, R = 2 Million 
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Figure 7.19 Upper Surface Velocity Profile, NACA 66 3 -018, 
X/C = .62, ACA = 0°, R = 2 Million 



84 



DISTANCE ADOVE SURFACE ETA/DELTA 



VELOCITY PROFILES 




Figure 7.20 



Upper Surface Velocity Profil 
X/C = .64, AOA =0°, R = 2 Million 
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Figure 7.21 Upper Surface Velocity Profile, NACA 66 3 -018, 
X/C = .66, AOA =0°, R = 2 Million 
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Figure 7.22 Upper Surface Velocity Profile, NACA 66 3 -018, 
X/C = .68, AOA = 0®, R = 2 Million 
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Figure 7.23 Upper Surface Velocity Profile, NACA 66 3 -018, 
X/C = .69, AOA = 0°, R = 2 Million 
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Figure 7.24 



Upper Surface Velocity Profile, NACA 663 - 018 , 
X/C = .70, AOA =0°, R = 2 Million 
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Figure 7.25 Upper Surface Velocity Profile, NACA 66 3 -018, 
X/C = .58, ACA = 2°, R = 2 Million 
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Figure 7.26 



Upper Surface Velocity Profile, NACA 66 3 -018, 
X/C = .60, AOA =2°, R = 2 Million 
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Figure 7.27 



Upper Surface Velocity Profile, NACA 66 3 -018, 
X/C = .62, AOA =2°, R = 2 Million 
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Figure 7.28 Upper Surface Velocity Profile, NACA 66 3 -018, 
X/C = .64, AOA =2°, R = 2 Million 
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Figure 7.29 



Upper Surface Velocity Profile, NACA 66 3 -C18, 
X/C = .66, AOA = 2°, R = 2 Million 
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Figure 7.30 Upper Surface Velocity Profile, NACA 66 3 -018, 
X/C = .68, AOA = 2°, R = 2 Million 
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Figure 7.31 



Upper Surface Velocity Profile, NACA 66 3 -018, 
X/C = .69, AOA =2°, R = 2 Million 
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In an effort to increase the region of separated flow, 
both the GY transition constant and the location of upper 
surface transition were adjusted. Table 7.1 shows 
theresults. With a constant of 1200 the transition location 
could not be moved aft. However, when the constant was 
lowered to 200 and below, the transition location, x/c = .69, 
found by Gault, could be used moving the transition inside 
the bubble. While a lowered GY tr . constant and increased x/c 
transition improved the bubble size, the separation length, 
x/c = .60 to .70, could not quite be met. The best result, 
bubble = .6391 - .7098 (x/c), was obtained with a GY tr . of 40, 
and an upper surface transition location, XTRU, input of .69 
(x/c). Whether 40 is a suitable value for other foils has 
not been determined. 

Twenty iterations were used on all computer runs for this 
airfoil. To make sure that 20 iterations were sufficient for 
accurate results, Figure 7.32 was obtained. The lift 
coefficient was plotted for each iteration, and as can be 
seen the even iterations produced very minimal changes past 
12. Even between 19 and 20 the change in lift was only 5 x 
10“*. Therefore, 20 iterations were considered sufficient 
for all computer runs. 

B. NACA 0010 (MODIFIED) 

Similar to the NACA 66 3 ~018, computer results of the NACA 
0010 (Modified) airfoil section were compared to the test 
results of Gault [Ref. 8]. The tests were also conducted in 
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TABLE 7.1 EFFECT OF GYTR AND XTRU ON THE LENGTH OF THE 
SEPARATION BUBBLE 



GYTR 


XTRU (X/C) SEPARATION (X/C) 


1200 


.620792 (1) 


6572 - 


.6929 




.65 


(2) 






.69 


(2) 




400 


.69 


(2) 




300 


.69 


(2) 




200 


.69 


.6391 - 


.7596 


120 


.639164 (1) 


No Separation 




.65 


.6572 - 


.6751 




.69 


.6391 - 


.7434 


80 


. 69 


.6391 - 


.7268 


60 


.69 


.6391 - 


.7268 


40 


.69 


.6391 - 


.7098 


30 


.69 


.6572 - 


.7098 


20 


.69 


.6572 - 


.7098 


10 


. 69 


.6572 - 


.7098 


NACA 663-018 








Reynolds Number = 


2,000,000 






Angle of Attack = 


0 Degrees 






Region of Separation, Experimental 


(Gault ) 


= .60 - .70 (X/C) 


GYTR = Empirical 


Transition Constant 




XTRU = Upper Surface Begin of Transition, 


Input or Computer 


Derived 








Lower Surface Begin of Transition 


= Computer Derived, Each 



Case 

1 ) Computer Derived 

2) Breakdown in Simulation 
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Figure 7.32 Lift Coefficient Versus Iterations, NACA 

66 3 -018, AOA = 0, R = 2 Million, Transition 
Constant = 1200 
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the 7-by-10-foot wind tunnel at NASA Ames and two-dimensional 
flow was simulated. 

The 0010 (modified) airfoil shown in Figure 7.33, unlike 
the 663-OI8, was not computer tested for section lift 
coefficients as Abbott and Doenhoff [Ref. 9] had no 
comparable airfoil . 

Leading edge, upper surface, laminar to turbulent, 
transition locations were observed, though, and the results 
are shown in Figures 7.34 and 7.35 for Reynolds numbers of 
two and six million, respectively. The Cebeci Code curves 
represent the beginning of transition. 

Full chord pressure distributions for angles of attack of 
zero and three degrees and Reynolds numbers of three and 
eight million are plotted in Figures 7.36 through 7.39. 

Leading edge pressure distributions for angles of attack 
of four, eight and twelve degrees, and Reynolds numbers of 
two and six million are shown in Figures 7.40 through 7.45. 
Of particular interest are the "lump" disparities in Figures 
7.42 through 7.45. A possible explanation for the computer 
program deletions of the lumps is a failure to predict 
leading edge bubbles. 

C. NACA 4412 

Computer results of the NACA 4412 airfoil section were 
compared to the test results of Hastings and Williams [Ref. 
10], which were performed in the 13-by 9-foot low speed wind 
tunnel of the Royal Aircraft Establishment at Bedford. The 
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Figure 7.34 Upper Surface Laminar to Turbulent 
Transition, Leading Edge, NACA 0010 
(Modified) , R = 2 Million 
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Figure 7.35 Upper Surface Laminar to Turbulent 
Transition, Leading Edge, NACA 0010 
(Modified) , R = 6 Million 
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Figure 7.36 Upper Surface Pressure Distribution, NACA 
0010 (Modified), AOA = 0°, R = 3 Million 
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Figure 7.37 



Upper Surface Pressure Distribution, NACA 
0010 (Modified), AOA = 0°, R = 8 Million 
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Figure 7.38 Upper Surface Pressure Distribution, NACA 
0010 (Modified), AOA = 3°, R = 3 Million 
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igure 7.39 Upper Surface Pressure Distribution, NACA 
0010 (Modified), AOA = 3°, R = 8 Million 
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Figure 7.40 Leading Edge Upper Surface Pressure 

Distribution, NACA 0010 (Modified) , AOA = 4°, 
R = 2 Million 
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Figure 7.41 Leading Edge Upper Surface Pressure 

Distribution, NACA 0010 (Modified), AOA = 4°, 
R = 6 Million 
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Figure 7.42 Leading Edge Upper Surface Pressure 

Distribution, NACA 0010 (Modified), AOA = 8°, 
R = 2 Million 
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Figure 7.43 Leading Edge Upper Surface Pressure 

Distribution, NACA 0010 (Modified), AOA = 8°, 
R = 6 Million 
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Figure 7.44 Leading Edge Upper Surface Pressure 
Distribution, NACA 0010 (Modified), 
AOA = 12°, R = 2 Million 
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Figure 7.45 Leading Edge Upper Surface Pressure 
Distribution, NACA 0010 (Modified) , 
AOA = 12° , R = 6 Million 
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one meter chord model spanned the full width, 13 feet, to 
simulate two-dimensional flow. 

Mounted at its quarter-chord point, the model was 
extensively instrumented with static pressure orifices. 
Boundary layer and wake measurements were made at mid-span 
where the 88 pressure orifices were located. 

The main emphasis in the experiment was on defining the 
upper surface boundary layer through separation and into the 
wake. Laser anemometry was used to measure the average 
velocities and Reynolds stresses were measured by hot-wire 
anemometry . 

The 4412 shown in Figure 7.46 was initially computer 
tested with the Cebeci code for momentum thickness. In 
Figure 7.47 the upper and lower surface laminar to turbulent 
transitions were computer derived, x/c t r- upper and lower 
surfaces = 0.00625. In Figure 7.48 the upper and lower 
surface transitions were input as x/c C r- upper surface = 0.01, 
and x/Ctr lower surface = 0.11. These values were as close 
as could be input to 0.014 and 0.110 respectively, for the 
downstream ends of the transition trips used in the 
experiment. The differences in transition locations seemed 
to make no difference in computer results. The momentum 
thicknesses still did not agree very well with Hastings and 
Williams' experimental results. 

Figure 7.49 compares lift coefficients from the Cebeci 
code, Hastings and Williams, and Abbott and Doenhoff [Ref. 
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Figure 7.46 NACA 4412 
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Figure 7.47 Upper Surface Momentum Thickness, NACA 4412, 
AOA = 12.49°, R = 4.17 Million, Computer 
Derived Transitions 
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Figure 7.48 Upper Surface Momentum Thickness, NACA 4412, 
AOA = 12.49°, R = 4.17 Million, Computer 
Derived Input Transitions 
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Figure 7.49 Lift Coefficients, NACA 4412 
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10]. The dashed line for the Cebeci code, Reynolds number 
equal to 4.17 million and computer derived transitions, 
should lie between the curves for Reynolds number equal to 
three and six million from Abbott and Doenhoff. However, it 
lies above the six million curve which further indicates an 
insufficient boundary layer development. With the input 
transition the code prediction for Reynolds number equal to 
4.17 million and angle of attack equal to 12.49 degrees, was 
even higher. Of interesting note though, is that the 
Hastings and Williams prediction, Reynolds number equals 4.17 
million and angle of attack equals 12.49 degrees, lies below 
the Abbott and Doenhoff values, possibly indicating an error 
on their part. 

Figures 7.50 through 7.56 show the velocity profiles from 
x/c = .66 to the trailing edge. In all cases the Reynolds 
number was 4.17 million, the angle of attack was 12.49 
degrees, and the upper and lower transitions were .01 and 
.11, respectively. U/U e indicates the fraction of the 
velocity at the boundary layer edge, and /Delta is the 
fraction of boundary layer thickness, where Delta, 6, is 
defined as the layer thickness where the velocity is 99% of 
the edge velocity. 

As these figures indicate, as well as Figure 7.57, a 
Cebeci code velocity profile summation, the code does predict 
separation, but not the extent indicated by Hastings and 
Williams. If the lift coefficient curves, Figure 7.58, can 



119 



00 10 - 



UPPER SURFACE VELOCITY PROFILE 




FT A /nFTTA 

NACA 4412 R=4.17 MIL A0A=12.49 

Figure 7.50 Upper Surface Velocity Profile, NACA 4412, 
X/C = .66, AOA = 12.49°, R = 4.17 Million 
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Figure 7.51 Upper Surface Velocity Profile, NACA 4412, 
X/C = .74, AOA = 12.49°, R = 4.17 Million 
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igure 7.52 Upper Surface Velocity Profile, NACA 4412, 
X/C = .78, AOA = 12.49°, R = 4.17 Million 
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Figure 7.53 Upper Surface Velocity Profile, NACA 4412, 
X/C = .85, AOA = 12.49°, R = 4.17 Million' 
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Figure 7.54 Upper Surface Velocity Profile, NACA 4412, 
X/C = .92, AOA = 12.49°, R = 4.17 Million 
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Figure 7.55 Upper Surface Velocity Profile, NACA 4412, 
X/C = .95, AOA = 12. 49°, R = 4.17 Million 
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Figure 7.56 Upper Surface Velocity Profile, NACA 4412, 
X/C = .997, AOA = 12.49°, R = 4.17 Million 
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Figure 7.58 Lift Coefficient, NACA 4412, 20 Iterations, 
R = 1.523 Million 
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be a reference, then it appears that the Cebeci code predicts 
an underdeveloped flow, too little separation, and Hastings 
and Williams show an overdeveloped flow, too much separation, 
for the given conditions. 

To insure that the Cebeci Code was run correctly, certain 
results by Cebeci, Clark, Chang, Halsey and Lee [Ref. 1] were 
attempted to be duplicated. Figure 7.58 compares two curves 
from Figure 14, [Ref. 1], curves labeled interactive theory 
and interactive theory with a modified transition, with a 
curve obtained using the Cebeci Code with 20 iterations and 
computer derived transitions. Interestingly, the first and 
third curves of Figure 7.58, interactive theory and Cebeci 
Code, respectively, should be the same, but the two clearly 
are not above nine degrees angle of attack. Even more 
interestingly, the Cebeci Code lift curve in Figure 7.59 
after only 10 iterations does match the interactive theory 
curve . 

Figure 7.60 clearly shows the importance of using enough 
iterations to obtain a reasonably accurate solution. 

Finally, Figure 7.61 shows a very good match between the 
pressure coefficients for set conditions of Figure 16, Cebeci 
et al [Ref. 1], and the Cebeci Code, 20 iterations. 

D. FX 63-137 

Computer results of the Wortmann FX 63-137 airfoil were 
compared to the test results of Brendel and Mueller [Ref. 
11], which were conducted in the University of Notre Dame 
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Figure 7.59 Lift Coefficient, NACA 4412, 20 Iterations, 
R = 1.523 Million 
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Figure 7.60 Lift Coefficient Versus Iterations, 

NACA 4412, AOA = 12°, R = 1.523 Million 
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Figure 7.51 



Variation of Pressure Coefficient, 

NACA 4412, AOA = 12°, R = 1.523 Million 
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.61m x .61m wind tunnel. Two cast epoxy resin airfoil models 
with chords of .305m and spans of ,4m were mounted in the 
center of the test section. Pressure was recorded on one 
model with 96 pressure taps connected through two scanivalves 
to an electronic manometer. Boundary layer velocity 
measurements were obtained on the other model using a 
constant temperature anemometer with a five pm diameter, 
single-sensor, hot-wire, boundary layer probe. 

Using the Cebeci Code the FX 63-137 airfoil shown in 
Figure 7.62 was initially tested for section lift 
coefficients with a transition constant of 1200. Reynolds 
numbers of .28, .5 and .7 million were used, and the results, 
shown in Figures 7.63 and 7.64, were compared to those of 
Althaus and Wortmann [Ref. 12]. Interestingly, the Cebeci 
Code predicted low values for Reynolds numbers of .28 and .5 
million, but for .7 million the lift coefficients were nearly 
identical to Althaus and Wortmann up to an angle of attack of 
10 degrees. 

As the purpose of Brendel and Mueller was to make 
boundary layer measurements at low Reynolds numbers, a 
computer comparison was unsuccessfully attempted for steady 
flow at a Reynolds number of 100,000 and an angle of attack 
of 7 degrees. With 20 iterations the Cebeci Code failed. 

To understand why the code calculations ceased for this 
case, other computer runs were attempted for the same 
Reynolds number and angle of attack, but with fewer 
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Figure 7.62 Wortmann FX 63-137 
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Figure 7.63 Lift Coefficient, FX 63-137, R = .28 Million, 
and .7 Million 



135 



0.2 0.4 0.6 0.8 



LIFT COEFFICIENT 



o 




Figure 7.64 Lift Coefficient, FX 63-137, R = .5 Million 
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iterations. Figures 7.65 and 7.66 show the upper surface 
displacement and momentum thicknesses for steady flow and 
iterations of 2, 4, 6,8 and 10. As can be seen in both 
figures flow calculations matched very well with experimental 
data up to approximately x/C = .55. After that point stall 
occurred and calculations ceased with more than 10 
iterations. Brendel and Mueller experimentally derived 
separation to begin at x/C = .34, but reattachment was shown 
to occur at x/c = .60. Unfortunately the Cebeci Code could 
not predict reattachment for the prescribed conditions. 
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Figure 7.65 Upper Surface Displacement Thickness, 
FX 63-137, AOA = 7°, R = 100,000 
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Figure 7.66 



Upper Surface Momentum Thickness, FX 63-137, 
AOA = 7°, R = 100,000 
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VIII. CONCLUSIONS AND RECOMMENDATIONS 



Cebeci's viscous/inviscid interaction program was applied 
to the analysis of steady, two dimensional, incompressible 
flow past four airfoils, the NACA 663 -OI 8 , 0010 (modified), 
4412 and the Wortmann FX 63-137. Detailed comparisons with 
the available experimental results show that for attached 
flows the essential features are correctly modelled, but that 
significant discrepancies are found in regions of flow 
separation. These discrepancies are possibly caused by the 
empirical transition modelling used in the present code. 
Future efforts therefore should be directed to the 
incorporation of transition calculations which permit the 
prediction of transition within a separation bubble, such as 
the application of the e^-method proposed by Cebeci [Ref. 
13]. 
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